function [clusterCenters normalizationFactors priorProb obsProb transProb] = extractFeaturesAndFormClusters( masterFolder, validFolders, models, handles )
    %First extract features from all images in the masterFolder subfolders
    %Feature extraction is done in the way specified by the model
    %Cluster this data set into M central vector (again M specified by the
    %model
    
    %Number of models
    nmodels = size(models, 2);
    %Alphabet size - current is 26 only capital letters
    alphabetSize = 26;
    %All training data gathered distinctly for each model
    allObservations = cell(nmodels, alphabetSize);
    rawObservations = cell(nmodels,1);
    for i = 1 : nmodels
        rawObservations{i} = [];
    end
    %Number of states in each model
    nmodelStates = zeros(nmodels,1);
    %Features - equals to 2 atm
    nfeatures = 2;
    
    for i = 1 : size(validFolders)
        %Construct the folder path from master folder path and valid folder
        %name
        folder_path = strcat(masterFolder,'\', validFolders(i,:));
        %Get the valid images in that folder
        validImages = getValidImages(folder_path);
        nimages = size(validImages,1);
        
        %Get which letter we are training
        trainingLetter = validFolders(i,1);
        %Get the number of states for this letter (each model)
        nmodelStates = getLetterStateNumber( models, trainingLetter);
        %Allocate storage for observations indexable
        allObservations{:,i} = zeros(nimages,nmodelStates,nfeatures);
            
        %Set the plotting structure, currently it is 3 image per rows
        plotRows = ceil(nimages/3);
        %Read and process images iteratively
        for j = 1 : nimages
            %Read the image
            img = imread(strcat(folder_path,'\', validImages(j,:)));
            %Threshold to a binary image
            img = threshold2Bin(img);
            %Get the letter boundary in the image
            [rmin rmax cmin cmax] = bboxletter(img,1);
            %Image cropped to letter area
            letterImg = img(rmin:rmax, cmin:cmax);
            %Plot the image
            lh = subplot(plotRows,3, j,'Parent',handles.pnlAllTrainingImagesInFolder);
            imshow(letterImg, 'Parent', lh);
            observations = letterObservations(letterImg, models, handles, 1, trainingLetter);
            
            for st = 1 : nmodels
                allObservations{st,i}(j,:,:) = observations{st};
                rawObservations{st} = [rawObservations{st};observations{st}];
            end
            
            %Check the execution method and wait for user click if necessary
            isPausedBetweenImages = get(handles.rbImages, 'Value');
            if isPausedBetweenImages
                %Wait for user to press the button
                waitfor(handles.btnTrain, 'UserData', 1);
                %Set the value back to 0 -corresponds to not pressed yet-
                set(handles.btnTrain, 'UserData', 0);
            end
        end

        %Check the execution method and wait for user click if necessary
        isPausedBetweenLetters = get(handles.rbLetters, 'Value');
        if isPausedBetweenLetters
            %Wait for user to press the button
            waitfor(handles.btnTrain, 'UserData', 1);
            %Set the value back to 0 -corresponds to not pressed yet-
            set(handles.btnTrain, 'UserData', 0);
        end
        
        %Then clear the axes
        delete(get(handles.pnlAllTrainingImagesInFolder, 'Children'));
        
    end
    
    clusteredObservations = cell(nmodels, alphabetSize);
    
    %Normalize the columns of raw data so each feature has similar effect
    %in the clustering
    for i = 1 : nmodels
        %Sum all observations for each feature distinctly
        featureTotals = sum(rawObservations{i});
        %Make the first feature reference. Its mean will be the mean of
        %other features also
        normalizationFactors{i} = featureTotals/featureTotals(1);
        %Normalize the data
        normalizedObservations{i} = rawObservations{i}./repmat(normalizationFactors{i},[size(rawObservations{i},1) 1]);
        %Cluster normalized data into specified number of clusters (M is
        %the number of observations that is going to be used in HMM)
        %Empty action, start and replicate are for safety. All of them
        %partly assures correct clustering
        [~,clusterCenters{i}] = kmeans(normalizedObservations{i}, models{i}.M, 'emptyaction', 'singleton', 'replicates', 5, 'start', 'uniform');
        
        
        %Now change our previous raw observations to cluster indexes
        for j = 1 : alphabetSize
            
            %Get which letter we are training
            trainingLetter = validFolders(j,1);
            %Get the number of states for this letter (each model)
            nmodelState = getLetterStateNumber( models, trainingLetter, i);
            
            %Construct the folder path from master folder path and valid folder
            %name
            folder_path = strcat(masterFolder,'\', validFolders(j,:));
            %Get the valid images in that folder
            nimages = size(getValidImages(folder_path),1);
            
            %Allocate storage
            clusteredObservations{i,j} = zeros(nimages,nmodelState);
            
            %Convert raw data to clusters
            for k = 1 : nimages
                clusterIndices = getCluster(allObservations{i,j}(k,:,:), clusterCenters{i});
                clusteredObservations{i,j}(k,:) = clusterIndices';
            end
            
            %Set prior probabilities of states
            %Trivial - only first state is nonzero
            priorProb{i,j} = zeros(nmodelState, 1);
            priorProb{i,j}(1) = 1;
            
            %Set initial guesses for transition probabilities
            transProb{i,j} = setInitialTransitionProbabilities(nmodelState);
            
            %Count the number of occurences of observations (cluster
            %indices) for each state
            obsProb{i,j} = zeros(nmodelState, models{i}.M);
            for l = 1 : nmodelState
                %Start from one for each observation frequency to compromise 
                %for the small training set we have
                 tempFreq = histc(clusteredObservations{i,j}(:,l), (1:models{i}.M)) + 1 ;
                 %Convert from frequency to probability
                 obsProb{i,j}(:,l) = tempFreq / sum(tempFreq);
            end
        end 
    end
end

